**************************************************************
*** Programs to calculate model fit and direct/indirect effects in SAR.
***
*** Created: 5-12-15
*** Modified: 9-5-19
***
**************************************************************

mata: mata clear



**************************************************************
****************** PROGRAMS **********************************
**************************************************************
version 10.1
mata:
mata set matastrict on 
void calceffsar(real matrix draws, real matrix eye, real matrix dub, real scalar sims)
{
	real matrix pd, effect, o1, o2, o3, dub2, dub3
	real scalar i, teffect, teffect_se, deffect, deffect_se, ieffect, ieffect_se, feffect, feffect_se, heffect, heffect_se, deffect0, deffect0_se, teffect1, teffect1_se, deffect1, deffect1_se, ieffect1, ieffect1_se, teffect2, teffect2_se, deffect2, deffect2_se, ieffect2, ieffect2_se, teffect3, teffect3_se, deffect3, deffect3_se, ieffect3, ieffect3_se
	
	dub2 = dub*dub
	dub3 = dub*dub*dub
	
	effect = J(sims,14,.)
	for(i=1; i<=sims; i++) {
		pd = luinv(eye-draws[i,1]*dub) * draws[i,2]
		
		effect[i,1]=1/rows(dub)*sum(pd)
		effect[i,2]=1/rows(dub)*trace(pd)
		effect[i,3]=1/rows(dub)*(sum(pd) - trace(pd))
		
		o1 = draws[i,1]*dub*draws[i,2]
		
		effect[i,4]=1/rows(dub)*sum(o1)
		effect[i,5]=1/rows(dub)*trace(o1)
		effect[i,6]=1/rows(dub)*(sum(o1) - trace(o1))

		o2 = draws[i,1]^2*dub2*draws[i,2]
		
		effect[i,7]=1/rows(dub)*sum(o2)
		effect[i,8]=1/rows(dub)*trace(o2)
		effect[i,9]=1/rows(dub)*(sum(o2) - trace(o2))
		
		o3 = draws[i,1]^3*dub3*draws[i,2]
		
		effect[i,10]=1/rows(dub)*sum(o3)
		effect[i,11]=1/rows(dub)*trace(o3)
		effect[i,12]=1/rows(dub)*(sum(o3) - trace(o3))	
		
		effect[i,13]=1/rows(dub)*(trace(pd) - sum(eye*draws[i,2]))
		effect[i,14]=1/rows(dub)*(sum(pd) - sum(o1) - sum(eye*draws[i,2]))
	}
	
	teffect = mean(effect[,1])
	teffect_se = sqrt(variance(effect[,1]))
	deffect = mean(effect[,2])
	deffect_se = sqrt(variance(effect[,2]))
	ieffect = mean(effect[,3])
	ieffect_se = sqrt(variance(effect[,3]))

	deffect0 = mean(draws[,2])
	deffect0_se = sqrt(variance(draws[,2]))
	
	teffect1 = mean(effect[,4])
	teffect1_se = sqrt(variance(effect[,4]))
	deffect1 = mean(effect[,5])
	deffect1_se = sqrt(variance(effect[,5]))
	ieffect1 = mean(effect[,6])
	ieffect1_se = sqrt(variance(effect[,6]))

	teffect2 = mean(effect[,7])
	teffect2_se = sqrt(variance(effect[,7]))
	deffect2 = mean(effect[,8])
	deffect2_se = sqrt(variance(effect[,8]))
	ieffect2 = mean(effect[,9])
	ieffect2_se = sqrt(variance(effect[,9]))

	teffect3 = mean(effect[,10])
	teffect3_se = sqrt(variance(effect[,10]))
	deffect3 = mean(effect[,11])
	deffect3_se = sqrt(variance(effect[,11]))
	ieffect3 = mean(effect[,12])
	ieffect3_se = sqrt(variance(effect[,12]))

	feffect = mean(effect[,13])
	feffect_se = sqrt(variance(effect[,13]))
	heffect = mean(effect[,14])
	heffect_se = sqrt(variance(effect[,14]))
	
	st_matrix("r(mateffect)", effect)

	st_numscalar("r(teffect)", teffect)
	st_numscalar("r(teffect_se)", teffect_se)
	st_numscalar("r(deffect)", deffect)
	st_numscalar("r(deffect_se)", deffect_se)
	st_numscalar("r(ieffect)", ieffect)
	st_numscalar("r(ieffect_se)", ieffect_se)	

	st_numscalar("r(deffect0)", deffect0)
	st_numscalar("r(deffect0_se)", deffect0_se)
	
	st_numscalar("r(teffect1)", teffect1)
	st_numscalar("r(teffect1_se)", teffect1_se)
	st_numscalar("r(deffect1)", deffect1)
	st_numscalar("r(deffect1_se)", deffect1_se)
	st_numscalar("r(ieffect1)", ieffect1)
	st_numscalar("r(ieffect1_se)", ieffect1_se)	

	st_numscalar("r(teffect2)", teffect2)
	st_numscalar("r(teffect2_se)", teffect2_se)
	st_numscalar("r(deffect2)", deffect2)
	st_numscalar("r(deffect2_se)", deffect2_se)
	st_numscalar("r(ieffect2)", ieffect2)
	st_numscalar("r(ieffect2_se)", ieffect2_se)	
	
	st_numscalar("r(teffect3)", teffect3)
	st_numscalar("r(teffect3_se)", teffect3_se)
	st_numscalar("r(deffect3)", deffect3)
	st_numscalar("r(deffect3_se)", deffect3_se)
	st_numscalar("r(ieffect3)", ieffect3)
	st_numscalar("r(ieffect3_se)", ieffect3_se)
	
	st_numscalar("r(feffect)", feffect)
	st_numscalar("r(feffect_se)", feffect_se)
	st_numscalar("r(heffect)", heffect)
	st_numscalar("r(heffect_se)", heffect_se)
}
end


version 10.1
mata:
mata set matastrict on 
void calceffslx(real matrix draws, real matrix dub, real scalar sims)
{
	real matrix o, effect
	real scalar i, teffect, teffect_se, deffect, deffect_se, ieffect, ieffect_se

	effect = J(sims,3,.)
	for(i=1; i<=sims; i++) {
		o = draws[i,1]*dub
		
		effect[i,1]=1/rows(dub)*sum(o)
		effect[i,2]=1/rows(dub)*trace(o)
		effect[i,3]=1/rows(dub)*(sum(o) - trace(o))
		
	}

	teffect = mean(effect[,1])
	teffect_se = sqrt(variance(effect[,1]))
	deffect = mean(effect[,2])
	deffect_se = sqrt(variance(effect[,2]))
	ieffect = mean(effect[,3])
	ieffect_se = sqrt(variance(effect[,3]))

	st_matrix("r(mateffect)", effect)

	st_numscalar("r(teffect)", teffect)
	st_numscalar("r(teffect_se)", teffect_se)
	st_numscalar("r(deffect)", deffect)
	st_numscalar("r(deffect_se)", deffect_se)
	st_numscalar("r(ieffect)", ieffect)
	st_numscalar("r(ieffect_se)", ieffect_se)
}
end

version 10.1
mata:
mata set matastrict on 
void calceffslx1(real matrix draws, real matrix dub1, real scalar sims)
{
	real matrix o1, effect
	real scalar i, teffect, teffect_se, deffect, deffect_se, ieffect, ieffect_se, deffect_o0, deffect_se_o0, teffect_o1, teffect_se_o1, deffect_o1, deffect_se_o1, ieffect_o1, ieffect_se_o1

	effect = J(sims,3,.)
	for(i=1; i<=sims; i++) {
		o1 = draws[i,2]*dub1
		
		effect[i,1]=1/rows(dub1)*sum(o1)
		effect[i,2]=1/rows(dub1)*trace(o1)
		effect[i,3]=1/rows(dub1)*(sum(o1) - trace(o1))
	}

	teffect_o1 = mean(effect[,1])
	teffect_se_o1 = sqrt(variance(effect[,1]))
	deffect_o1 = mean(effect[,2])
	deffect_se_o1 = sqrt(variance(effect[,2]))
	ieffect_o1 = mean(effect[,3])
	ieffect_se_o1 = sqrt(variance(effect[,3]))
	
	deffect_o0 = mean(draws[,1])
	deffect_se_o0 = sqrt(variance(draws[,1]))
	
	teffect = mean(draws[,1] + effect[,1])
	teffect_se = sqrt(variance(draws[,1] + effect[,1]))
	deffect = mean(draws[,1] + effect[,2])
	deffect_se = sqrt(variance(draws[,1] + effect[,2]))
	ieffect = mean(effect[,3])
	ieffect_se = sqrt(variance(effect[,3]))
	
	st_matrix("r(mateffect)", effect)

	st_numscalar("r(teffect_o1)", teffect_o1)
	st_numscalar("r(teffect_se_o1)", teffect_se_o1)
	st_numscalar("r(deffect_o1)", deffect_o1)
	st_numscalar("r(deffect_se_o1)", deffect_se_o1)
	st_numscalar("r(ieffect_o1)", ieffect_o1)
	st_numscalar("r(ieffect_se_o1)", ieffect_se_o1)
	
	st_numscalar("r(deffect_o0)", deffect_o0)
	st_numscalar("r(deffect_se_o0)", deffect_se_o0)
	
	st_numscalar("r(teffect)", teffect)
	st_numscalar("r(teffect_se)", teffect_se)
	st_numscalar("r(deffect)", deffect)
	st_numscalar("r(deffect_se)", deffect_se)
	st_numscalar("r(ieffect)", ieffect)
	st_numscalar("r(ieffect_se)", ieffect_se)
}
end

version 10.1
mata:
mata set matastrict on 
void calceffslx2(real matrix draws, real matrix dub1, real matrix dub2, real scalar sims)
{
	real matrix o1, o2, effect
	real scalar i, teffect, teffect_se, deffect, deffect_se, ieffect, ieffect_se, feffect, feffect_se, heffect, heffect_se, deffect_o0, deffect_se_o0, teffect_o1, teffect_se_o1, deffect_o1, deffect_se_o1, ieffect_o1, ieffect_se_o1, teffect_o2, teffect_se_o2, deffect_o2, deffect_se_o2, ieffect_o2, ieffect_se_o2

	effect = J(sims,8,.)
	for(i=1; i<=sims; i++) {
		o1 = draws[i,2]*dub1
		o2 = draws[i,3]*dub2
		
		effect[i,1]=1/rows(dub1)*sum(o1)
		effect[i,2]=1/rows(dub1)*trace(o1)
		effect[i,3]=1/rows(dub1)*(sum(o1) - trace(o1))
		
		effect[i,4]=1/rows(dub2)*sum(o2)
		effect[i,5]=1/rows(dub2)*trace(o2)
		effect[i,6]=1/rows(dub2)*(sum(o2) - trace(o2))
		
		effect[i,7]=1/rows(dub1)*trace(o2)
		effect[i,8]=1/rows(dub1)*sum(o2)
	}

	deffect_o0 = mean(draws[,1])
	deffect_se_o0 = sqrt(variance(draws[,1]))
	
	teffect_o1 = mean(effect[,1])
	teffect_se_o1 = sqrt(variance(effect[,1]))
	deffect_o1 = mean(effect[,2])
	deffect_se_o1 = sqrt(variance(effect[,2]))
	ieffect_o1 = mean(effect[,3])
	ieffect_se_o1 = sqrt(variance(effect[,3]))

	teffect_o2 = mean(effect[,4])
	teffect_se_o2 = sqrt(variance(effect[,4]))
	deffect_o2 = mean(effect[,5])
	deffect_se_o2 = sqrt(variance(effect[,5]))
	ieffect_o2 = mean(effect[,6])
	ieffect_se_o2 = sqrt(variance(effect[,6]))
	
	feffect = mean(effect[,7])
	feffect_se = sqrt(variance(effect[,7]))
	heffect = mean(effect[,8])
	heffect_se = sqrt(variance(effect[,8]))
	
	teffect = mean(draws[,1] + effect[,1] + effect[,4])
	teffect_se = sqrt(variance(draws[,1] + effect[,1] + effect[,4]))
	deffect = mean(draws[,1] + effect[,2] + effect[,5])
	deffect_se = sqrt(variance(draws[,1] + effect[,2] + effect[,5]))
	ieffect = mean(effect[,3] + effect[,6])
	ieffect_se = sqrt(variance(effect[,3] + effect[,6]))
	
	st_matrix("r(mateffect)", effect)

	st_numscalar("r(teffect_o1)", teffect_o1)
	st_numscalar("r(teffect_se_o1)", teffect_se_o1)
	st_numscalar("r(deffect_o1)", deffect_o1)
	st_numscalar("r(deffect_se_o1)", deffect_se_o1)
	st_numscalar("r(ieffect_o1)", ieffect_o1)
	st_numscalar("r(ieffect_se_o1)", ieffect_se_o1)
	
	st_numscalar("r(teffect_o2)", teffect_o2)
	st_numscalar("r(teffect_se_o2)", teffect_se_o2)
	st_numscalar("r(deffect_o2)", deffect_o2)
	st_numscalar("r(deffect_se_o2)", deffect_se_o2)
	st_numscalar("r(ieffect_o2)", ieffect_o2)
	st_numscalar("r(ieffect_se_o2)", ieffect_se_o2)
	
	st_numscalar("r(deffect_o0)", deffect_o0)
	st_numscalar("r(deffect_se_o0)", deffect_se_o0)
	
	st_numscalar("r(feffect)", feffect)
	st_numscalar("r(feffect_se)", feffect_se)
	st_numscalar("r(heffect)", heffect)
	st_numscalar("r(heffect_se)", heffect_se)
	
	st_numscalar("r(teffect)", teffect)
	st_numscalar("r(teffect_se)", teffect_se)
	st_numscalar("r(deffect)", deffect)
	st_numscalar("r(deffect_se)", deffect_se)
	st_numscalar("r(ieffect)", ieffect)
	st_numscalar("r(ieffect_se)", ieffect_se)
}
end

version 10.1
mata:
mata set matastrict on 
void calceffslx3(real matrix draws, real matrix dub1, real matrix dub2, real matrix dub3, real scalar sims)
{
	real matrix o1, o2, o3, effect
	real scalar i, teffect, teffect_se, deffect, deffect_se, ieffect, ieffect_se, feffect, feffect_se, heffect, heffect_se, deffect_o0, deffect_se_o0, teffect_o1, teffect_se_o1, deffect_o1, deffect_se_o1, ieffect_o1, ieffect_se_o1, teffect_o2, teffect_se_o2, deffect_o2, deffect_se_o2, ieffect_o2, ieffect_se_o2, teffect_o3, teffect_se_o3, deffect_o3, deffect_se_o3, ieffect_o3, ieffect_se_o3

	effect = J(sims,11,.)
	for(i=1; i<=sims; i++) {
		o1 = draws[i,2]*dub1
		o2 = draws[i,3]*dub2
		o3 = draws[i,4]*dub3
		
		effect[i,1]=1/rows(dub1)*sum(o1)
		effect[i,2]=1/rows(dub1)*trace(o1)
		effect[i,3]=1/rows(dub1)*(sum(o1) - trace(o1))
		
		effect[i,4]=1/rows(dub2)*sum(o2)
		effect[i,5]=1/rows(dub2)*trace(o2)
		effect[i,6]=1/rows(dub2)*(sum(o2) - trace(o2))

		effect[i,7]=1/rows(dub3)*sum(o3)
		effect[i,8]=1/rows(dub3)*trace(o3)
		effect[i,9]=1/rows(dub3)*(sum(o3) - trace(o3))	
		
		effect[i,10]=1/rows(dub2)*(trace(o2) + trace(o3))
		effect[i,11]=1/rows(dub2)*(sum(o2) + sum(o3))
	}

	deffect_o0 = mean(draws[,1])
	deffect_se_o0 = sqrt(variance(draws[,1]))
	
	teffect_o1 = mean(effect[,1])
	teffect_se_o1 = sqrt(variance(effect[,1]))
	deffect_o1 = mean(effect[,2])
	deffect_se_o1 = sqrt(variance(effect[,2]))
	ieffect_o1 = mean(effect[,3])
	ieffect_se_o1 = sqrt(variance(effect[,3]))

	teffect_o2 = mean(effect[,4])
	teffect_se_o2 = sqrt(variance(effect[,4]))
	deffect_o2 = mean(effect[,5])
	deffect_se_o2 = sqrt(variance(effect[,5]))
	ieffect_o2 = mean(effect[,6])
	ieffect_se_o2 = sqrt(variance(effect[,6]))
	
	teffect_o3 = mean(effect[,7])
	teffect_se_o3 = sqrt(variance(effect[,7]))
	deffect_o3 = mean(effect[,8])
	deffect_se_o3 = sqrt(variance(effect[,8]))
	ieffect_o3 = mean(effect[,9])
	ieffect_se_o3 = sqrt(variance(effect[,9]))
	
	feffect = mean(effect[,10])
	feffect_se = sqrt(variance(effect[,10]))
	heffect = mean(effect[,11])
	heffect_se = sqrt(variance(effect[,11]))
	
	teffect = mean(draws[,1] + effect[,1] + effect[,4] + effect[,7])
	teffect_se = sqrt(variance(draws[,1] + effect[,1] + effect[,4] + effect[,7]))
	deffect = mean(draws[,1] + effect[,2] + effect[,5] + effect[,8])
	deffect_se = sqrt(variance(draws[,1] + effect[,2] + effect[,5] + effect[,8]))
	ieffect = mean(effect[,3] + effect[,6] + effect[,9])
	ieffect_se = sqrt(variance(effect[,3] + effect[,6] + effect[,9]))
	
	st_matrix("r(mateffect)", effect)

	st_numscalar("r(deffect_o0)", deffect_o0)
	st_numscalar("r(deffect_se_o0)", deffect_se_o0)
	
	st_numscalar("r(teffect_o1)", teffect_o1)
	st_numscalar("r(teffect_se_o1)", teffect_se_o1)
	st_numscalar("r(deffect_o1)", deffect_o1)
	st_numscalar("r(deffect_se_o1)", deffect_se_o1)
	st_numscalar("r(ieffect_o1)", ieffect_o1)
	st_numscalar("r(ieffect_se_o1)", ieffect_se_o1)
	
	st_numscalar("r(teffect_o2)", teffect_o2)
	st_numscalar("r(teffect_se_o2)", teffect_se_o2)
	st_numscalar("r(deffect_o2)", deffect_o2)
	st_numscalar("r(deffect_se_o2)", deffect_se_o2)
	st_numscalar("r(ieffect_o2)", ieffect_o2)
	st_numscalar("r(ieffect_se_o2)", ieffect_se_o2)

	st_numscalar("r(teffect_o3)", teffect_o3)
	st_numscalar("r(teffect_se_o3)", teffect_se_o3)
	st_numscalar("r(deffect_o3)", deffect_o3)
	st_numscalar("r(deffect_se_o3)", deffect_se_o3)
	st_numscalar("r(ieffect_o3)", ieffect_o3)
	st_numscalar("r(ieffect_se_o3)", ieffect_se_o3)
	
	st_numscalar("r(feffect)", feffect)
	st_numscalar("r(feffect_se)", feffect_se)
	st_numscalar("r(heffect)", heffect)
	st_numscalar("r(heffect_se)", heffect_se)
	
	st_numscalar("r(teffect)", teffect)
	st_numscalar("r(teffect_se)", teffect_se)
	st_numscalar("r(deffect)", deffect)
	st_numscalar("r(deffect_se)", deffect_se)
	st_numscalar("r(ieffect)", ieffect)
	st_numscalar("r(ieffect_se)", ieffect_se)
}
end


version 10.1
mata:
mata set matastrict on 
void calceff2(real matrix draws, real matrix eye, real matrix dub, real scalar sims)
{
	real matrix pd1, pd2, effect
	real scalar i, teffect1, teffect1_se, teffect2, teffect2_se, deffect1, deffect1_se, ieffect1, ieffect1_se, deffect2, deffect2_se, ieffect2, ieffect2_se
	
	effect = J(sims,6,.)
	for(i=1; i<=sims; i++) {
		/*** Beta_1 ***/
		pd1 = luinv(eye-draws[i,1]*dub) * draws[i,2]
		
		effect[i,1]=1/rows(dub)*sum(pd1)
		effect[i,2]=1/rows(dub)*trace(pd1)
/*		_diag(pd1,0) */
		effect[i,3]=1/rows(dub)*(sum(pd1) - trace(pd1))
		
		/*** Beta_2 ***/
		pd2 = luinv(eye-draws[i,1]*dub) * draws[i,3]
		
		effect[i,4]=1/rows(dub)*sum(pd2)
		effect[i,5]=1/rows(dub)*trace(pd2)
/*		_diag(pd2,0)	*/
		effect[i,6]=1/rows(dub)*(sum(pd2) - trace(pd2))
	}

	teffect1 = mean(effect[,1])
	teffect1_se = sqrt(variance(effect[,1]))
	deffect1 = mean(effect[,2])
	deffect1_se = sqrt(variance(effect[,2]))
	ieffect1 = mean(effect[,3])
	ieffect1_se = sqrt(variance(effect[,3]))

	teffect2 = mean(effect[,4])
	teffect2_se = sqrt(variance(effect[,4]))
	deffect2 = mean(effect[,5])
	deffect2_se = sqrt(variance(effect[,5]))
	ieffect2 = mean(effect[,6])
	ieffect2_se = sqrt(variance(effect[,6]))
	
	st_matrix("r(mateffect)", effect)

	st_numscalar("r(teffect_x1)", teffect1)
	st_numscalar("r(teffect_se_x1)", teffect1_se)
	st_numscalar("r(deffect_x1)", deffect1)
	st_numscalar("r(deffect_se_x1)", deffect1_se)
	st_numscalar("r(ieffect_x1)", ieffect1)
	st_numscalar("r(ieffect_se_x1)", ieffect1_se)

	st_numscalar("r(teffect_x2)", teffect2)
	st_numscalar("r(teffect_se_x2)", teffect2_se)	
	st_numscalar("r(deffect_x2)", deffect2)
	st_numscalar("r(deffect_se_x2)", deffect2_se)
	st_numscalar("r(ieffect_x2)", ieffect2)
	st_numscalar("r(ieffect_se_x2)", ieffect2_se)	
}
end

capture program drop spreg_fit
program define spreg_fit, rclass
syntax
	preserve
		local Nobs = e(N)
		local dv = e(depvar)
		keep in 1/`Nobs'
		
		* Calculate the measures of fit.
		sum `dv', meanonly
		scalar y_mean = r(mean)	
	
		* Total (spatial + non-spatial) model fit
		predict y_hat in 1/`Nobs', rform
		egen double mss = sum((y_hat - y_mean)^2)
		egen rss = sum((`dv' - y_hat)^2)
	
		* Non-spatial model fit (XB)
		predict y_hat_xb in 1/`Nobs', xb
		egen double mss_xb = sum((y_hat_xb - y_mean)^2)
		egen rss_xb = sum((`dv' - y_hat_xb)^2)
	
		egen double tss = sum((`dv' - y_mean)^2)
	
		keep in 1

		* Total (spatial + non-spatial) model fit
		gen double r2 = mss/tss
		gen double rmse = sqrt((rss)/e(N))
		local p = e(k) - e(df_m)
		gen double ar2 = 1-(1-r2)*((e(N)-1)/(e(N)-`p'-1))

		sum r2, meanonly
		return scalar r2 = round(r(mean), .0001)
		sum ar2, meanonly
		return scalar ar2 = round(r(mean), .0001)
		sum rmse, meanonly
		return scalar rmse = round(r(mean), .0001)
	
		* Non-spatial model fit (XB)
		gen double r2_xb = mss_xb/tss	
		gen double rmse_xb = sqrt((rss_xb)/e(N))
		local p_xb = e(k) - e(df_m)
		gen double ar2_xb = 1-(1-r2_xb)*((e(N)-1)/(e(N)-`p_xb'-1))
	
		sum r2_xb, meanonly
		return scalar r2_xb = round(r(mean), .0001)
		sum ar2_xb, meanonly
		return scalar ar2_xb = round(r(mean), .0001)
		sum rmse_xb, meanonly
		return scalar rmse_xb = round(r(mean), .0001)

		* AIC and BIC
		estat ic
		mat S = r(S)
		return scalar aic = round(S[1,5], .001)
		return scalar bic = round(S[1,6], .001)
	
	restore
end

capture program drop sr_fit
program define sr_fit, rclass
syntax
	preserve

	local dv = e(depvar)
	
	predict y_hat, xb
	
	* Calculate the measures of fit.
	qui sum `dv', meanonly
	scalar y_mean = r(mean)	
	gen mu_hat = (`dv' - y_hat)^2
	gen num = (y_hat - y_mean)^2
	gen denom = (`dv' - y_mean)^2
	collapse (sum) ess = num tss = denom rss = mu_hat
	gen r2 = ess/tss
	gen rmse = sqrt(rss/e(N))
	local p = e(k) - e(df_m)
	gen adj_r2 = 1-(1-r2)*((e(N)-1)/(e(N)-`p'-1))
	qui sum r2, meanonly
	return scalar r2 = r(mean)
	qui sum adj_r2, meanonly
	return scalar ar2 = r(mean)
	qui sum rmse, meanonly
	return scalar rmse = r(mean)
	
	* AIC and BIC
	local ll = e(ll)
	local p = e(rank0) + e(k_eq_model)
	local N = e(N)
	local aic = (-2 * `ll') + (2 * `p')
	return scalar aic = `aic'
	local bic = (-2 * `ll') + (`p' * ln(`N'))
	return scalar bic = `bic'

	capture drop cons y_hat1 mu_hat num denom r2 rmse adj_r2 
	restore
end

**************************************************************
**************************************************************
**************************************************************
